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Abstract. - The response of turbulent flow to time-modulated forcing is studied by direct 
numerical simulations of the Navier-Stokes equations. The large-scale forcing is modulated via 
periodic energy input variations at frequency u). The response is maximal for frequencies in 
the range of the inverse of the large eddy turnover time, confirming the mean-field predictions 
of von der Heydt, Grossmann and Lohse (Phys. Rev. E 67, 046308 (2003)). In accordance 
with the theory the response maximum shows only a small dependence on the Reynolds number 
and is also quite insensitive to the particular flow-quantity that is monitored, e.g., kinetic energy, 
dissipation-rate, or Taylor- Reynolds number. At sufficiently high frequencies the amplitude 
of the kinetic energy response decreases as 1/uj. For frequencies beyond the range of maximal 
response, a significant change in phase-shift relative to the time-modulated forcing is observed. 



Introduction. Recently, response maxima in time modulated turbulence have been 
predicted within a mean field theory of turbulence [1] . Subsequently, such response maxima 
were found [2] in numerical simulations of simplified dynamical turbulence models such as 
the GOY model [3-5] or the reduced wave vector approximation (REWA) [6-9]. However, 
these response maxima computed in [2] were not pronounced at all, due to the approximate 
treatment of the small scales in either of these approaches. Indications of response maxima 
resulting from time-modulated forcing have subsequently also been seen in experiment [10]. 
The experimental observations were done by introducing a time-dependent swirl to fluid in a 
closed container and monitoring the energy-dissipation-rate. The selected set-up did not allow 
to identify possible flow-structuring under resonance conditions, nor to conclusively distinguish 
such resonance phenomena from flow-organization associated with the size of the container. 
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The purpose of this paper is to complement these theoretical, numerical, and experimental 
observations by direct numerical simulations (DNS) of turbulence, subject to time-modulated 
large-scale forcing. In a turbulent flow whose large-scale forcing is periodically modulated in 
time, all typical flow-properties develop a complex time-dependence. However, averaging such 
turbulent time-dependence, conditioned on the phase of the periodic modulation, yields a clear 
and much simpler periodic pattern [2] . The dependence of the conditionally averaged response 
on the frequency of the modulation may be quantified by monitoring changes in flow-properties 
such as total energy, dissipation-rate, or Taylor- Reynolds number. In case of a fast modulation 
with a frequency w > w^, where lul is the inverse large eddy turnover time, only a modest 
effect on the flow is expected, or none at all. Likewise, if lo <C ojl the modulation is quasi- 
stationary and the flow may be expected to closely resemble the corresponding unmodulated 
case. In between these extremes a more pronounced response may develop, which is the subject 
of this investigation. 

The DNS approach allows to investigate in detail the response of turbulent flow-properties 
to periodic modulation of the forcing. In particular, we present an extensive parameter- 
study involving a large range of modulation frequencies for two different Reynolds numbers, 
and establish response maxima in a variety of flow-properties. The response is found to be 
significantly increased at modulation frequencies on the order of the inverse of the eddy- 
turnover time. Near resonance, the 'activity' of the turbulent flow is found to be considerably 
higher than in the unmodulated case. At high frequencies uj the amplitude of the modulation- 
specific response of the kinetic energy is found to uniformly decrease to zero as a; -1 . This type 
of external control of turbulence may offer new opportunities with relevance to technological 
applications. 

The organization of this paper is as follows. We first introduce the computational flow- 
model in more detail. Subsequently, an overview of the ensemble averaging procedure and data 
extraction is given. Then the main result, the response of various flow properties to time- 
modulated forcing, is presented. The paper ends with a summary and conclusions. 

Computational flow-model. The full Navier-Stokes equations for incompressible flow 
are numerically solved in a periodic flow-domain with a pseudo-spectral code. In spectral 
space, the Navier-Stokes equations read 



u a (k,t) = M al 3y(k) ^2 ^(P^Wfa)*) + F a (k,t), (1) 

p+q=k 



with M a £ 7 (k) = y % \ kpD a ^(k) + k 7 D a 0(k)J , with D a p(k) = 5 a p — k a kp/\k\ 2 . Here, v is 

the kinematic viscosity, u a (k,t) is the Fourier-coefficient of the velocity field at wave vector 
k and time t and F a is the time-modulated forcing. 

First, we recall that traditional agitation of the large-scale structures in a turbulent flow 
may be achieved by introducing a forcing term restricted to wave vectors with |k| < kp, i.e., 
identifying a forcing-range through the upper- limit hp. Specifically, we force the turbulence 
similarly as in [7, 11], 

f^) = WV7^ ; |k|<fc ^ (2) 

jVf |u(k, t)\ 

where e w is the constant energy injection rate and Np = Np(kp) is the total number of forced 
modes. For convenience, the wave vectors are grouped in spherical shells with the n-th shell 
containing all modes such that (n — 1/2) < |k| < (n + 1/2). We applied large-scale forcing 
either in the first shell at n = 1 (i.e., kp = 3/2 which implies Np = 18, the case considered 
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in [2]) or in the first two shells (i.e., kp = 5/2 which implies Np — 80). The second step in 
specifying the forcing F a introduces the periodic time modulation 



where Ap is the amplitude of modulation and u> its angular frequency. The modulated forcing 
corresponds to a total energy input rate which oscillates around e w with amplitude Ap, 



The lengths and time scales of the numerical simulation are chosen by picking L = 1 for 
the box-size in physical space and e w = 0.15 for the energy injection rate. The Reynolds 
number is then determined by the dimcnsionless viscosity v. Choosing v = 1060.7 and 
v^ 1 = 4242.6 result in respective approximate Taylor-Reynolds numbers R\ = 50 and 
R\ = 100. We use these two cases as references denoted by R50 and i?ioo- The spatial 
resolution needed may be estimated by requiring k max r/ > 1 [12] with r\ the Kolmogorov dissi- 
pation scale and k max the highest wavenumber included in the spatial discretization. For R 50 
case a resolution of at least iV 3 = 64 3 computational points is required while for i?ioo a higher 
resolution of 192 3 points is necessary. The latter poses a strong computational challenge in 
view of the extensive ensemble averaging and large number of modulation frequencies. How- 
ever, it was found that many large scale quantities, such as the total energy, do not depend too 
sensitively on resolution. As an example, a resolution of 64 3 points corresponds to k max r\ w 0.4 
for the i?ioo case. Still, this resolution is quite adequate for studying the response of total 
energy This was verified by repeating the analysis at a selection of modulation frequencies 
with resolutions 128 3 and 192 3 . The predictions of quantities that rely more on small scales, 
such as the dissipation-rate, contain a higher numerical uncertainty for i?ioo case and 64 3 com- 
putational points, but still allow a clear interpretation of the main turbulence response. This 
was separately assessed using the higher resolution data at selected characteristic frequencies. 

The direct numerical simulation for the unmodulated case starts from an initial condition 
that was generated on the basis of the Pao spectrum [13]. We adopt exactly the same initial 
conditions as in [14] which allow a separate validation of the simulations. Explicit second 
order compact-storage Runge-Kutta time-stepping [15] with fully de-aliased pseudo-spectral 
discretization is used. The unmodulated turbulent flows provide the point of reference for 
investigating the effect of modulated forcing, to which we turn next. 

Averaging procedure and simulation setting. - In order to analyze the response to a time- 
modulated forcing, the precise extraction of the amplitude and phase of the conditionally 
averaged variations is a key issue. Two steps can be distinguished, i.e., the computation 
of the conditionally averaged signal itself and the subsequent determination of amplitude and 
phase characteristics of this signal, see Figure These steps are discussed and illustrated 
next. 

We adopt ensemble averaging to determine the conditionally averaged signal S(uj,t), where 
S(u>,t) is the total energy E(w,t), the Taylor-Reynolds number R\(uj,t) or the energy dis- 
sipation rate e(cu,t). Ensemble averaging requires a sufficiently large sample of statistically 
independent signals {Sj(u>,t)} to be generated. The computational approach is summarized 
in Figure ^a) and involves two main steps. Firstly, we compute the unmodulated flow and 
store N r realizations of the turbulent solution corresponding to t > 10. The latter condition 
allows transients related to the initial condition to become negligible. The time-separation 
between these snapshots is larger than two eddy-turnover times. Subsequently, each of these 
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Fig. 1 - Procedure of data generation to compute the conditional average (a) and the amplitudes 
Af(u>), Ae(u>), A £ {uj) and phase-shifts <bp(uf) = 0, $ £ (a>) of the forcing Tp{u},t) (dashed 

line), the energy QE(u,t) (labeled o), and the energy-dissipation-rate Q e {u),t) (labeled o) normalized 
by their respective means Tf = s w , Q E (u>), and Q e (w) obtained from simulations at the modulation 
frequency uj — 0.8n (b). 



N r realizations was taken as the initial condition for a simulation with time-modulated forc- 
ing at a particular frequency lo. This provides N r sample signals which need to be averaged 
to obtain the conditionally averaged signal S(u),t). Repeating this procedure for a range 
of frequencies yields the total response characteristics. Given the conditionally averaged re- 
sponse signal S(u),t), there are various ways in which amplitude and phase information can 
be extracted. In [10] the signal S(u},t) was first averaged over time to yield S(uj). Subse- 
quently, the normalized variation defined as Qg (ui,t) = S(u),t)/S(u>) was studied using the 
Fourier-transform (J-) in which time t is transformed into frequency /. Correspondingly, the 
power amplitude spectrum Qg(u,f) = T^Q^^jt) — l) can be obtained which assumes 

a maximum value As(uj) = max{|Qg (cj, f)\}\f=f A (u)i as denoted in Figure^b) for forcing 
Af(u>), total energy Ae(lo), and energy-dissipation-rate A e {uj). The maximum A$(uj) as the 
amplitude at dominant frequency can be used to quantify the response as function of the 
modulation frequency to. This approach is accurate if Fourier-transform is applied to an inte- 
ger number of modulation periods. The method used in [2] is based on a fitting procedure in 
which it is assumed that S(u>, t) tts S + As sin (wt + ^s). The dependence of the parameters 
{S, As, $s} on u> may be obtained from a least squares procedure. This evaluation method 
assumes that the conditionally averaged signal has the same frequency as the forcing. 

At modest ensemble-size N r it is beneficial to explicitly incorporate variations in the un- 
modulated reference signal to improve the data-evaluation. This motivates an alternative 
method in which we determine N r sample signals {Sj(u>, t)} corresponding to the modulated 
case, as well as N r unmodulated signals {sj(t)} that start from the same set of initial condi- 
tions. This allows to generate different 'normalized' signals such as Qg(u>, t) = J2j Sjl 12 j s j 

or Qg(u),t) = Sj/sj/N r . These normalized signals provide estimates that compensate to 
some degree for the relatively small number of samples or for an unknown mean component 
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but have the drawback that they cannot be applied in the context of a physical experiment. 
Additionally, we divided these signals by its means (time-averages) and removed the con- 
stant component corresponding to the zero-frequency response. Application of the Fourier- 
transform, = J~(Q$ /Qs — l) and Q$ = T (Q^ / Q$ — l), provides direct access to 
amplitude and phase information. Each of these methods identified above yields the same 
general impression of response maxima in time-modulated turbulence. Differences arise only 
on a more detailed level of the processed data but these do not obscure the interpretation 
of the main features of the response. Therefore we only present results extracted from the 

normalized signal Qs/Qs = QsV*9s^ m wna t follows, unless explicitly stated otherwise. The 
simulations were performed in the frequency range ir/5 < lo < 80ir with time-modulated 
forcing at an amplitude Ap = 1/5. Referring to Figure ^ for each of the N r unmodulated 
initial conditions, tit = 4 periods of the modulated forcing were simulated, i.e., each sample 
signal was computed for tltT time- units with modulation-period T = 2tt/uj. Since an explicit 
time-stepping method was adopted, the cases at low to add particularly to the total computa- 
tional cost. The number of realizations required in the ensemble was investigated separately. 
Results for several modulation frequencies were compared at N r = 10, 30 and 50; it was found 
that 30 independent samples provide adequate statistical convergence for our purposes. We 
stored N t = 40 points per modulation period and present results obtained by evaluating the 
last two recorded periods, i.e., 2T < t < AT. Comparison with results obtained by evaluating 
data on < t < AT yielded only minor differences. Finally, the phase <&s(w) between the 
forcing and the response can be computed from the Fourier-transformed data as well. At the 
dominant frequency /a of the transformed signal Qs(ui,f) = J-(Qs(u, ty/Qsi^? — 1)j the 
phase becomes $s(w) = arctan (Iir(Qs{uj, /a))/Rc(Qs(o;, /a)))- 




Fig. 2 - The response QE(u,t)/Q E (u>) for the R50 case recorded at different modulation frequencies u> 
is shown in (a) together with the modulation of the forcing TF(u>,t)/e m (dashed). The corresponding 
power spectra of the Fourier-transform as function of the transformed frequency / are collected in 
(b). Modulation frequencies lo/(2tt) = 0.1, 0.2, 0.3, 0.4, 0.5, 2.0 are included and labeled by o, 0, □, X, >, 
and *, respectively. 
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Modulated turbulence. - In Figure|2Ia) the conditionally averaged signal Qe{u, t)/Q E (u>) 
based on total energy is shown at a number of modulation frequencies. The conditionally 
averaged response has a clear oscillatory behavior. The Fourier-transform of the data from 
Figure Ela) is shown in Figure Elb) and displays a dominant maximum corresponding to 
the forcing frequency /a = u>/(2tt). This observation confirms that the least-squares fitting 
procedure adopted in [2] is justified. 
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Fig. 3 - (a) Amplitude of the energy response Ab(w) and (b) compensated energy response ujAe(u>) 
obtained for the R50 (labeled o) and the -R100 (labeled >) cases. Verification at selected frequencies 
for resolution 128 3 and R100 case (labeled *). Results for forcing in two first shells (ftp < 5/2) and 
the R50 case (labeled o). The inset in (a) shows the phase-shift <3>_e(o;) between the energy response 
and the forcing modulation. 

We now focus on the amplitude of the total energy response as function of the modula- 
tion frequency to. The amplitude Ae(oS) computed as maximum of the Fourier-transformed 
normalized signal for each modulation frequency is shown in Figure |3Ja). The maximum 
response appears at u m ax ~ 1-5, in accordance with the expectation [1,2] that it should be 
close to the inverse large eddy turnover time. In addition, the location of the maximum is 
not very sensitive to R\, reflecting that the response maximum is mainly associated with 
the large-scale features in the flow. At high modulation frequencies u> > u moI the decay of 
Ae is proportional to uj^ 1 , which becomes particularly visible in the compensated response 
loAe(u), Figure[2Jb). At very low modulation frequencies u> < u> max a plateau in Ae{uj) must 
of course develop [1,2], as the turbulence then completely follows the forcing. Our simulations 
do not achieve small enough u) to observe a pronounced plateau. 

The maximum of loAe{lu) is about 35% higher as compared to the value at high cu. This is 
as expected lower than predicted by the mean-field theory described in [1] as the fluctuations 
slightly smear out the mean-field maximum, but it is much more pronounced compared to 
results based on the GOY or REWA simulations [2]. The reason is that, although the ap- 
pearance of the response maxima is a large-scale effect, the correct resolution of the small- 
scales is important for a proper quantitative representation of the effect, because the small 
scale resolution affects the energy flux downscale. We also calculated the response curves for 
the Taylor-Reynolds number; the results are quite similar. 

The phase-difference between the forcing modulation and the conditionally averaged total 
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Fig. 4 - (a) Amplitude of the energy-dissipation-rate response A e (u>) and (b) compensated energy- 
dissipation-rate response uiA e (uj) obtained for the R$o case. The inset in (a) shows the phase-shift 
$ £ (cj) between the energy-dissipation response and the forcing modulation. 



energy response is shown in Figure OH a) as inset. We observe a strong variation in this 
phase-difference for modulation frequencies near the most responsive modulation frequency. 
It appears that the maximum response as shown in FigureODoccurs at a modulation frequency 
where also the variation in the phase-difference is largest. A strong phase shift was found 
similarly in windtunncl experiments in which a time-modulation is introduced via a periodic 
cycling of an upstream active grid. In these experiments the maximum response was found to 
shift to higher frequencies in case the characteristic length-scales of the forcing were reduced. 

Can such a dependence on the type of forcing also be observed in our numerical sim- 
ulations? To find out we force a higher wavenumber band of modes {kp < 5/2) instead 
of restricting us entirely to low wavenumber forcing. The result is seen in Figure OHb) indi- 
cated by diamonds. Indeed, for this type of forcing the response maximum is less pronounced. 
Further quantitative connections with physical experiments are currently being investigated. 

The energy-dissipation-rate in the system is a quantity that is accessible to direct physical 
experimentation. In Figure 0] we show the energy-dissipation-rate amplitude A e (uj). We 
notice that at high modulation frequency lu the amplitude approaches zero, consistent with the 
expectation that the modulation of the forcing is not effective in this range. More importantly, 
the encrgy-dissipation-rate amplitude displays a strong response maximum at the level of 85% 
compared to the amplitude of modulation. The total mean energy-dissipation T -1 e(u>, t)dt 
for each modulation frequency u> is almost constant. It differs from the energy input rate 
e w = 0.15 at the level of 1% for most of the frequencies, reaching the maximum difference of 
5% for the lowest simulated frequency, confirming good numerical convergence. 

Summary and Conclusions. - The direct numerical simulation of the response of tur- 
bulence to time-modulated forcing confirms the existence of a response maximum. The sim- 
ulation findings are in general agreement with predictions based on a mean-field theory [1]. 
The mean-field theory predicts the decrease of the response amplitude proportional to uj^ 1 as 
the modulation frequency is sufficiently large which was observed in the simulations as well. 
The response maxima in the total energy and the Taylor-Reynolds number occur at the forcing 
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frequencies of the order of the inverse large eddy turnover time scale. The phase-difference be- 
tween the modulation of the forcing and the conditionally averaged response displays a strong 
dependence on the modulation frequency as well. The modulation frequency at which the 
response maximum arises depends only weakly on the Reynolds number but shows a depen- 
dence on the scales included in the forcing as well as on the flow-property that is considered. 
In general, if the particular quantity of interest shows a stronger dependence on the smaller 
scales in a turbulent flow, then the response maximum arises at a somewhat higher frequency. 
These findings may be independently assessed in physical experiments, e.g., conducted in wind 
tunnels combined with the use of active grids cycled in a periodic sequence [16]. 

* * * 

Stimulating discussions with Willem van de Water (Eindhoven University of Technology) 
are gratefully acknowledged. This work is part of the research program 'Turbulence and its 
role in energy conversion processes' of the Foundation for Fundamental Research of Matter 
(FOM), in the Netherlands, which is financially supported by the Netherlands Organization 
for Scientific Research (NWO). The authors wish to thank SARA Computing and Networking 
Services in Amsterdam for providing the computational resources. 

REFERENCES 

[1] A. von der Heydt, S. Grossmann, and D. Lohse, Phys. Rev. E 67, 046308 (2003). 

[2] A. von der Heydt, S. Grossmann, and D. Lohse, Phys. Rev. E 68, 066302 (2003). 

[3] L. Biferale, E. Calzavarini, F. Toschi, and R. Tripiccione, Europhys. Lett. 64, 461 (2003). 

[4] T. Bohr, M. H. Jensen, G. Paladin, and A. Vulpiani, Dynamical Systems Approach to Turbulence 
(Cambridge University Press, Cambridge, 1998). 

[5] L. Kadanoff, D. Lohse, J. Wang, and R. Benzi, Phys. Fluids 7, 617 (1995). 

[6] J. Eggers and S. Grossmann, Phys. Fluids A 3, 1958 (1991). 

[7] S. Grossmann and D. Lohse, Z. Phys. B 89, 11 (1992). 

[8] S. Grossmann and D. Lohse, Phys. Fluids 6, 611 (1994). 

[9] S. Grossmann and D. Lohse, Phys. Rev. E 50, 2784 (1994). 

[10] O. Cadot, J. H. Titon, and D. Bonn, J. Fluid Mech. 485, 161 (2003). 

[11] S. Ghosal, T. S. Lund, P. Moin, and K. Akselvoll, J. Fluid Mech. 286, 229 (1995). 

[12] V. Eswaran and S. B. Pope, Computers Fluids 16, 257 (1988). 

[13] S. B. Pope, Turbulent Flow (Cambridge University Press, Cambridge, 2000). 

[14] J. Meyers, B. J. Geurts, and M. Baelmans, Phys. Fluids 15, 2740 (2003). 

[15] B. J. Geurts, Elements of direct and large-eddy simulation (R.T. Edwards, 2003). 

[16] W. van der Water, 2005, private communication. 



